Exotic disordered phases in the quantum Ji — J2 model on the honeycomb lattice 
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We study the ground-state phase diagram of the frustrated quantum Ji — J2 Heisenberg antifer- 
romagnet on the honeycomb lattice using a mean field approach in terms of the Schwinger boson 
representation of the spin operators. We present results for the ground-state energy, local magne- 
tization, energy gap and spin-spin correlations. The system shows magnetic long range order for 
< J2/J1 < 0.2075 (Neel) and 0.398 < J2/J1 < 0.5 (spiral). In the intermediate region, we find 
two magnetically disordered phases: a gapped spin liquid phase which shows short-range Neel cor- 
relations (0.2075 < J2/J1 < 0.3732), and a lattice nematic phase (0.3732 < J2/J1 < 0.398), which 
is magnetically disordered but breaks lattice rotational symmetry. The errors in the values of the 
phase boundaries which are implicit in the number of significant figures quoted, correspond purely 
to the error in the extrapolation of our finite-size results to the thermodynamic limit. 

PACS numbers: 75.10.Kt, 75.10.Jm, 75.50.Ee 
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I. INTRODUCTION 

The two-dimensional (2D) Heisenberg model on bipar- 
tite lattices has been intensively studied in the last years. 
In the unfrustrated case, the classical ground state is ob- 
tained when all the spins in one sublattice are pointing 
in a given direction whereas in the other sublattice the 
spins are pointing in the opposite direction. However, 
in the quantum case this state is not the real ground 
state, in fact this is not an eigenstate of the Hamilto- 
nian. The quantum ground state is exactly known in one 
dimension^, but no exact results for the two dimensional 
antiferromagnet are known, even for simple lattices as 
the square lattice. However, several experimental and 
numerical studies suggested that the ground state is in 
fact the spin SU(2) symmetry broken Neel type state. In 
contrast, when we include frustration in the system, for 
example by including second neighbor interactions, the 
ground state may become much more complicated. 

In the quantum case, the ground state energy is lower 
than the classical value, due to the quantum fluctuations. 
The effects of these fluctuations vary depending on the 
dimension, the spin quantum number, the presence of 
frustrating interactions and the coordination number of 
the lattice. One can ask what the quantum fluctuations 
are when the coordination number is changed. In two 
dimensions two paradigmatic examples of unfrustrated 
systems are the square lattice, with coordination number 
z = A, and the honeycomb lattice with z = 3. Previous 
results^i^ have shown that the staggered magnetization 
is smaller in the z = 3 case. This behavior is in accord 
with the tendency towards a less classical behavior for 
systems of lower coordination number. 

The inclusion of frustration in 2D quantum antiferro- 
magnets is expected to enhance the effect of quantum 
spin fluctuations and hence suppress magnetic order—. 
This idea has motivated many researchers to look for 
its realization^"—. A special scenario to check this is the 
frustrated Heisenberg model on the honeycomb lattice. 
Due to the small coordination number (z = 3) which is 




FIG. 1. (Color online) The honeycomb lattice 

with Ji and J2 couplings considered in this paper. 
The lattice sites with different colors belong to differ- 
ent sublattices. The primitive translation vectors of 
the direct lattice are [ei= (^3/2, 3/2) , 62= (73/2, -3/2)] . 
ai=(0,-l), a2=(V3/2, 1/2) and 33= (-%/3/2, 1/2) corre- 
spond to the nearest neighbor bonds. 



the lowest allowed in a 2D system, quantum fluctuations 
could be expected to be stronger than those on the square 
lattice and may destroy the antiferromagnetic order— "J^. 

The study of frustrated quantum magnets on the hon- 
eycomb lattice has also experimental motivational^"—. 
One of the most exciting experimental progresses is one 
kind of bismuth oxynitrate, Bi3Mn40i2(N03), which was 
obtained by Smirnova et ali^. In this compound the 
Mn"'+ ions form a S* = 3/2 honeycomb lattice with- 
out any distortion. The magnetic susceptibility data in- 
dicates two-dimensional magnetism. Despite the large 
AF Weiss constant of -257K, no long-range ordering 
was observed down to 0.4K, which suggests a nonmag- 
netic ground statei^"— . The substitution of Mn^+ in 
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FIG. 2. (Color online) Phase diagram as a function of the 
frustration J2/J1. a) Classical phase diagram, b) Quantum 
phase diagram corresponding to 5* = i obtained by means of 
SBMFT. 



Bi3Mn40i2(N03) by V^"*" may lead to the realization of 
the 5* = 1/2 Heisenberg model on the honeycomb lattice. 

The analysis of the honeycomb lattice from a more 
general point of view has gained lately a lot of interest 
both coming from graphene-related issues^l and from the 
possible spin-liquid phase found in the Hubbard model in 
such gcomctry22ri2^. Due to these reasons, recently there 
is huge theoretical interest in frustrated Heisenberg mod- 
els on the honeycomb lattice, in which frustration is in- 
corporated by second nearest neighbors coupling a^'^'^"' '— 
and maybe also third nearest neighbors couplings^^"— . 

Motivated by previous results, in this paper we study 
the spin- 1/2 Heisenberg model on the honeycomb lattice 
with first (Ji) and second (J2) neighbors couplings. Us- 
ing a Schwinger boson mean field theory (SBMFT) we 
find strong evidence for the existence of an intermediate 
disordered region where a spin gap opens and spin-spin 
correlations decay exponentially. This magnetically dis- 
ordered region quantitatively agrees well with recent nu- 
merical simulation result o'^^'^^i'^^'^^ . Another key finding 
of our work is the presence of two kinds of magnetically 
disordered phases in this region. One is a gapped spin liq- 
uid (GSL) ^^'^^ with short-range Neel correlations, main- 
taining the lattice translational and rotational symmetry. 
The other phase is a staggered dimer valence-bond crys- 
tal (VBC), which is also called lattice nematio^S. This 
phase breaks lattice rotational symmetry, but preserves 
lattice translational symmetry. 

The rest of the paper is arranged as follows. In Sec. II 
we introduce our model and give a quick overview of the 
final phase diagram. In Sec. HI the general formalism 
of the Schwinger boson mean-field approach is presented. 
In Sec. IV, using the solutions of mean field equations, 
we discuss the phase diagram, especially the magneti- 
cally disordered region. We close with a summary and 
discussion in Sec. V. 




FIG. 3. (Color online) Sketch of the staggered dimer VBC 
state which breaks the lattice rotational symmetry but pre- 
serves the lattice translational symmetry. 



II. MODEL AND OVERVIEW OF THE PHASE 
DIAGRAM 

The Ji — J2 Heisenberg model on the honeycomb lattice 
is given by 

= Jl ^ Sx ■ Sy + J2 ^ Sx ■ Sy, (1) 

(xy>i <xy)2 

where Sx is the spin operator on site x and (xy)„ in- 
dicates sum over the n-th neighbors (see Fig. [1]). In 
this paper we are interested in the antiferromagnetic case 
{Ji, J2 > 0), and we focus on the region J2/J1 € [0,0.5]. 

In the classical limit, 5 — >■ c», the model displays dif- 
ferent zero temperature phases^"—, see Fig. [DJa). For 
J2/J1 < 1/6, the system is Neel ordered, while for 
J2/J1 > 1/61 the system shows spiral phases. For 
the quantum case, aspects of this model have been 
explored previously in the literature by various ap- 
proaches, including spin wave theory^i^i^i^, non-linear 
(T-modcl approach^, mean field theoryiSiS^^ exact diago- 
nalization (ED)^i^i^, variational Monte Carlo (VMC) 
method^i^, series expansion (SE)^, pseudofermion 
functional rcnormalization group (PFFRG)^ and cou- 
pled cluster method (CCM)^. However, these works 
yielded conflicting physical scenarios. 

This model was studied by Mattsson et ali^ using 
SBMFT with a mean field decoupling that considers 
only antiferromagnetic correlations for nearest neighbors 
and ferromagnetic correlations for next nearest neigh- 
bors. This scheme can only correctly describe Neel order. 
More recently Wang2^ studied this model within SBMFT 
including antiferromagnetic correlations for both near- 
est and next nearest neighbors. Unfortunately, The au- 
thor did not give the phase diagram for different val- 
ues of J2I J\- Actually, for frustrated models we can 
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FIG. 4. (Color online) Gap in the boson dispersion extrap- 
olated to the thermodynamic limit as a function of the frus- 
tration J2/J1 corresponding to S — 1/2. The gapped region 
corresponds to two different magnetically disordered phases: 
one is GSL, the other is staggered dimer VBC. Inset: Z3 order 
parameter defined in Eq. (|32[) . The onset of the VBC phase 
is determined by the value of J2/J1 where |V'| is non-zero (red 
arrows) 



not generally exclude either ferromagnetic or antiferro- 
magnctic correlations^ and is important to use a mean 
field decomposition that allows to include ferromagnetic 
an antiferromagnetic correlations in equal footing. An- 
other point is that both of them assume the bond mean 
fields are independent of the directions of bonds. There- 
fore, these two schemes can not describe the phases in 
which the lattice rotational symmetry has broken. Here 
we study the Hamiltonian ^ in the strong quantum 
limit (5' = 1/2) using a rotationally invariant version of 
this technique, which has proven successful in incorpo- 
rating quantum fluctuations^^i^"— . 

Our main results are summarized in Fig. E^b). The 
magnetic phase diagram is divided into four regions 1^ 
At small values of the frustrating coupling J2/J1, the 
system presents a Neel-like ground state. By increasing 
the frustration, we find at J2/JI — 0.2075 a continuous 
transition to a gapped spin liquid phase. When the value 
of the frustrating coupling exceeds J2/JI — 0.3732, we 
find a continuous transition into a staggered dimer VBC 
(lattice nematic) with broken Z3 symmetry (See Fig. [3]), 
which transforms at J2/J1 — 0.398 into a spiral phase. 



III. SCHWINGER BOSON MEAN-FIELD 
APPROACH 

It is well known that the SBMFT provides a natural 
description for both magnetically ordered and disordered 



FIG. 5. (Color online) Local magnetization determined by 
Eq. (|31|l extrapolated to the thermodynamic limit as a func- 
tion of the frustration J2/J1. The shaded region corresponds 
to the magnetically disordered phases. Insets correspond to 
the regions where the magnetization for Neel (left) and Spiral 
(right) phases becomes zero. 



phases based on the picture of the resonating valence 
bond states^i^"— . As a merit, this method does not start 
from any magnetic long range order for the ground state 
(in contrast to spin wave theory), which should emerge 
naturally if the Schwinger bosons condense at some mo- 
mentum vector—. At this momentum vector, the lowest 
excitation spectrum of the Schwinger bosons should be 
gapless. On the other hand. If the Schwinger bosons are 
gapped, the phase is magnetically disordered. In the fol- 
lowing, we will present in detail the rotationally invariant 
version of SBMFT which was introduced by Ceccatto et 
al^"— and we use in the following sections. 

Consider the SU(2) Heisenberg Hamiltonian on a gen- 
eral lattice: 



(2) 



xyct;? 



where x and y are the positions of the unit cells and 
vectors Va are the positions of each atom within the unit 
cell. Jap (x — y) is the exchange interaction between the 
spins located in x -I- and y -I- r^. 

In what follows we assume that the classical order can 
be parameterized as 

^yL+Yc ^ sin (x) (3) 

Sl+Y^ = (4) 

5'^+r„ = S C0SV3q(x), (5) 

with (pa (x) = Q • X + 0q , where Q is the ordering vector 
and Oa arc the relative angles between the classical spins 
inside each unit cell. 

The spin operators Sx on site x are represented by two 
bosons 6x(T (ct =tj I) 



= i bt • CT br 



br = 



5rt 



(6) 
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where a = {c7x,<7y,<Jz) are the Pauli matrices. Eq. ([6]) is 
a faithful representation of the algebra SU(2) if wc take 
into account the following local constraint 

25 = S1^6xt + 6l^i. (7) 
The exchange term can be expressed as 

Sx+r„ • Sy+r^=: bIij (x, y)Bai3 (x, y) : - Aj^^(x, y)ia;3 (x, y) , 

(8) 

where Aq,_^(x, y) and i3Q,_^(x, y) arc SU(2) invariants de- 
fined as 

i„,^(x,y) = i^a6(;^)6(«, (9) 

(T 

B^A^,y)^^Y.^bi^f^yl (10) 

(7 

with a l- The double dots (: O :) indicate the normal 
ordering of operator O. This decoupling is particularly 
useful in the study of magnetic systems near disordered 
phases, because it allows to treat antiferromagnetism and 
ferromagnetism in equal footing^"—. On the other hand, 
this scheme has been tested to obtain quantitatively quite 
accurate results which show excellent agreements with 

p^ p38. 54-56 

To construct a mean field Hamiltonian we perform the 
following Hartree-Fock decoupling 

(Sx+rc • Sy+rp)MF = [B*p(x-y)Bc,fj{^,y) 

- - y)Aap{x,y) + H.c] 

- ((Sx+r„ ■ Sy+r^)i\fF), (11) 
I 



where we have defined 

^:/3(x-y) = (il^(x,y)) (12) 
B:Mx-y) = (Bl^(x,y)) (13) 

and ( ) denotes the expectation value in the ground state 
at T = 0. It is convenient to change variables to R = 
x — y, and eliminating x in the sums wc obtain 

RyQ/3 I CT 

- {\B^AR)f-\AcA^)n}. 

The mean field Hamiltonian is quadratic in the boson 
operators and can be diagonalized. It is convenient to 
transform the operators to momentum space 

where Nc is the number of unit cells. After some algebra 
and using the symmetry properties: 

A^p(R) = -Ap^i-R) (15) 

we obtain the following form for the Hamiltonian 



E ■^-'^ W [\BM'^)\' - , (16) 

Ra/3 



where 



B 



R 

(k)=iE>^-/'W^"MI^)e"*'^''^"°"'"'^ 

R 



A 



(17) 
(18) 
(19) 



Now, we impose the constraint ([7]) in average over each 
sublatticc a by means of Lagrange multipliers A*-"^ 



H 



MF 



Hmf + H\ 



(20) 



with 



xa \ cr / 



(21) 



Using the symmetries (fTS]) we can see that both kinds 
of bosons (t, i) give the same contribution to the Hamil- 
tonian. Then, we can perform the sum over a to obtain 
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Hm. = \ E {(7f,(k) + A(")5.,)6tf + (7f,(-k) + A(")^„,)S1(36L«, - a (7„^,(k)Stf ^fe^g + 7«^,(k)6(f 6^^^ J } 



ka/; 
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Ro 



(22) 



It is convenient to introduce the Nambu spinor (k) 
(^Id^^, b_ki) where 

b_k, = (SL(-\6L(-),...,S!.(-^)) (23) 

and ric is the number of atoms in the unit ceh. Now, we 
can rewrite the Hamihonian into a compact form: 

i^A/f =E bt(k).i?(k).b(k) (24) 

k 

-(25+l)iV,EA(")-(i/MF), 

OL 

where the 2 ric x 2 dynamical matrix i^(k) is given by 



^' \ 7^Mk) 7f/3(k) + A(")5., 

To diagonahze the Hamihonian (|24p we need to per- 
form a para-unitary transformation of the matrix D(k) 
which preserves the bosonic commutation relations^. 
We can diagonahze the Hamiltonian by defining the new 
operators k ~ F h, where the matrix F satisfy 



(Ft)--3.(F)-=r3, r3^( Y -L )■ ^''^ 



With this transformation, the Hamihonian reads 

Hmf = E 4 • E(k) • at - (25 + l)iV, ^ A^") - (Hmf), 

k a 

(26) 

where 

E (k) = diag(cji (k) , . . . , (k) , (k) , . . . , (k) ) . (27) 

In terms of the original bosonic operators, the mean 
field parameters are 



^kt "-ki/ 

Hbitb'S)} (28) 

k 

„ g-.k(R+r„-r,)^5t(g5M^^| (29) 

and the constraint in the number of bosons can be written 
in the momentum space as 



-■ik(R,+rn — 



t(-)SW) + (6t(") 



-k4. 



>} 



r 



where Nc is the total number of unit cells and S is the 
spin strength. The mean field equations ([28]) and ([29|l 
must be solved in a self-consistent way together with the 
constraints (pO]) on the number of bosons. 

Finding numerical solutions involves finding the roots 
of the coupled nonlinear equations for the parameters A 
and B, plus the additional constraints to determine the 
values of the Lagrange multipliers A'"-*. We perform the 
calculations for finite but very large lattices and finally 
we extrapolate the results to the thermodynamic limit. 

We solve numerically for different values of the frus- 
tration parameter J2/J1 and with the values obtained 
for the MF parameters and the Lagrange multipliers we 
compute the energy and the new values for the MF pa- 
rameters. We repeat this self-consistent procedure until 
the energy and the MF parameters converge. After reach- 
ing convergence we can compute all physical quantities 
like the energy, the excitation gap, the spin-spin correla- 
tion and the local magnetization. During the calculation, 
it is convenient to fix the energy scale by setting the value 
of the nearest-neighbor coupling Ji = 1. 



IV. RESULTS 

In Fig. [4l we show the boson dispersion relation gap 
extrapolated to the thermodynamic limit as a function 
of the frustration (J2/J1). In the gapped region, the 
absence of Bose condensation indicates that the ground 
state is magnetically disordered. This result agrees well 
with recent ED^, VMO^e and CCMS>^ studies. In the 
gapless region, the excitation spectrum is zero at a given 
wave vector k* = Q/2, where the Boson condensation 
occurs. This is characteristic of the magnetically ordered 
phases. The structure of these phases can be understood 
through the spin-spin correlation function (SSCF) and 
the excitation spectrum. Some typical examples for dif- 
ferent phases will be shown later. 

To pin down the precise phase boundaries between the 
magnetically ordered and disordered phases, we intro- 
duce the local magnetization M(Q) as an order parame- 
ter, which is obtained from the long distance behavior of 
the spin-spin correlation function (SSCF)^!^: 



lim 



M2(Q)cos[Q.(x-y)]. (31) 



= 2SNc 



(30) 



In Fig. [51 we show the local magnetization for J2/J1 G 
[0,0.5]. For J2/J1 = 0, the local magnetization is 
AI{Q) =0.24176, which is in excellent agreement with 
the second order spin wave calculation result of 0.241 81^ 
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FIG. 6. (Color online) Ground-state energy per unit cell 
extrapolated to the thermodynamic limit as a function of the 
frustration J2/J1. The regions of the four different phases are 
indicated using the same colors that are used in Fig. [2l 

This value is significantly reduced by quantum fiuctua- 
tions compared with the classical value 0.5. The quantum 
Monte Carlo (QMC) result^ is 0.2677(6), which is con- 
siderably larger than ours. For the unfrustrated case, all 
the mean field approaches are quite inaccurate compared 
with much more controlled techniques like QMC. The 
difference in the M(Q) values of about 10%, provides, 
in the absence of any other quantitative evidence for the 
accuracy of the method as applied to this model, an indi- 
cation of the accuracy of the method and of all the results 
quoted that depend on the order parameters, including 
the phase boundaries. However, the mean field approach 
is still very useful to study gapped phases in frustrated 
systems. On one hand it is well known that for frus- 
trated systems QMC presents the famous sign problem. 
On the other hand, the study of quantities like energy 
gap requires the study of big sizes clusters and the use 
of exact diagonalization for small size clusters makes it 
very difficult to extrapolate the results. 

As J2I J\ increases, the local magnetization decreases. 
It vanishes continuously at J2I J\ — 0.2075, as shown in 
Fig. [5j^ This value is in excellent agreement with recent 
numerical results, such as 0.2 by Mezzacapo et alS^ us- 
ing VMC with an entangled-plaquette variational ansatz, 
as well as 0.207 ± 0.003 by Bishop et al^ using CCM. 
The shift of Neel boundary compared with the classical 
estimate 1/6 is due to quantum fluctuations which pre- 
fer to coUincar Neel rather than spiral phases in some 
casesi^ In this region, the SSCF is antiferromagnetic in 
all directions, and the Boson condensation happens at 
the r point of the first Brillouin zone: k* = (0, 0), which 
corresponds to the ordering vector Q = (0,0). As J2/J1 
decreases from 0.5, the local magnetization M(Q) de- 
creases. It vanishes continuously at J2/J1 — 0.398, as 
shown in Fig. This value is also in good agreement 
with recent numerical results, such as 0.4 by Mezzacapo 



et al^, as weU as 0.385 ± 0.010 by Bishop et alM. In 
this region, the SSCF shows different properties in differ- 
ent directions, however, it exhibits long range order in all 
directions. The gapless points of the excitation spectrum 
move continuously inside the first Brillouin zone as J2/ Ji 
changes. This results correspond to a spiral phase. In the 
classical version (S — 00) of the model (See Fig. [^IJa)), 
for J2/J1 > 1/6 there remains a line-type degeneracy in 
which the spiral wave number is not determined uniquely 
and is allowed on a ring in the Brillouin zone4^i^ Our 
results suggest that the classical degeneracy is lifted in 
the quantum version, where some spiral wave vectors are 
favored by quantum fluctuations from the manifold of 
classically degenerate spiral wave vectors. This spiral or- 
der by disorder selection was already seen by using a spin 
wave approach by Mulder et al. ^ and we have recovered 
this selection with a different approach. 

The most interesting part of the phase diagram is the 
intermediate region which has no classical counterpart. 
In this region, the nonmagnetic ground state retains 
SU(2) spin rotational symmetry and the lattice trans- 
lational symmetry. However, it may break the direc- 
tional symmetry of the lattice. Following Mulder et al^ 
we introduce the Z3 directional symmetry breaking order 
parameter \'ip\ where 

V = (Sa (r) • Ss (r)) + uj (S^ (r) • Ss (r + ei)) 

H-u;2 (r) • Sb (r - 62)) . (32) 

Here A, B correspond to the two different sublattices, 
r denotes the unit cell position, and w = exp(j27r/3). 
Equivalently, Okumura et al^ deflne 1113 = eiai-|-£2a2 + 
£333, where (11 = 1, 2, 3) are bond energies correspond- 
ing to the three nearest neighbor bonds = 1,2,3). 
It is trivial to see = |m3|. This order parameter 
is zero when the spin correlations along the three di- 
rections are equal. We flnd that \^\ keeps zero when 
•h/'h ^ 0.3732; it becomes non-zero continuously at 
J2/ Jl ~ 0.3732 as shown in Fig. Therefore, in the 
region 0.2075 < J2/J1 < 0.3732, the ground state pre- 
serves the Z3 lattice rotational symmetry. The SSCF 
shows short range antiferromagnetic correlations in all 
directions, and the minimum of the excitation spectrum 
remains pinned at the F point. Namely, the system re- 
mains to be a GSL. The appearance of the GSL agrees 
with recent two different VMC studiesi^i^ In the region 
0.3732 < J2/J1 < 0.398, the Z3 lattice rotational symme- 
try has broken. We find that the values of the mean fields 
A and B: A{B)i,.^ = ^(-S)a3 7^ A{B)e,^ \ the bond ener- 
gies have the same property: £2 = £3 7^ £i- Therefore, 
the system should belong to the staggered dimer VBC 
(lattice ncmatic). To further analyze this region, one 
need to calculate the dimer-dimer correlations. However, 
it is out of the scope of the present paper. The existence 
of the staggered dimer VBC is in agreement with a recent 
ED study,— a bond operator mean field study,— and a 
VMC study,^ 

The errors in the values of the phase boundaries that 
arc implicit here in the number of significant figures 
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FIG. 7. (Color online) SSCF for a system of size iV = 2 x 
50 X 50 in the zigzag direction corresponding to the four differ- 
ent phases: (a)J2/Ji = 0.18 (Neel), (b)J2/Ji = 0.36 (GSL), 
(c) J2/J1 = 0.38 (staggered dimer VBC), and (d) J2/J1 = 0.48 
(spiral) . 



FIG. 8. (Color online) SSCF for a system of size iV = 
2 X 50 X 50 in the armchair direction corresponding to the 
four different phases: (a)J2/Ji = 0.18 (Neel), (b)J2/Ji = 
0.36 (GSL), (c)J2/Ji = 0.38 (staggered dimer VBC), and 
(d)J2/Ji = 0.48 (spiral). 



quoted, correspond purely to the error in the extrapola- 
tion of our finite-size results to the thermodynamic limit. 
In no way are they intended to represent the essentially 
unknown errors implicit in the mean- field approach, e.g., 
the 10% difference in M (Q) compared with the QMC re- 
sult in the unfrustrated limit. All the transition values 
presented in this paper correspond to mean field estima- 
tions. In order to improve these values, it is necessary 
to study in detail the phase transitions beyond the mean 
field level, which is out of the scope of the present paper. 

In Fig. [6] we show the results for the ground state en- 
ergy per unit cell extrapolated to the thermodynamic 
limit. For the unfrustrated case (J2 = 0), Egs/Nc=- 
1.09779, which is in excellent agreement with the second 
order spin wave calculation result of — 1.0978i^ Com- 
pared with published QMC results by Reger et al^: 
-1.0890(9), and more recently by Low^: -1.08909(39), 
it has appreciable difference, as our previous discussion 
of the difference in the M(Q) values. Since energy esti- 
mates always have an intrinsic quadratic error, compared 
to an intrinsic linear error for other properties, even small 
errors in the energy can be of significance. The shape of 



the energy curve also supports that the three quantum 
phase transitions are continuous. 

In the following we show several typical examples for 
the four different phases. The SSCF along zigzag and 
armchair directions for a system of 5000 sites is shown 
in Fig. Hand Fig.[5]for J2/J1 = 0.18 (Neel), 0.36 (GSL), 
0.38 (staggered dimer VBC) and 0.48 (spiral). The corre- 
sponding lowest excitation spectrum is shown in Fig. [9l 
Although it is a finite size system, we can still see the 
corresponding properties for the four different phases as 
we have presented above. For J2/J1 = 0.18, the SSCF 
in both of the zigzag and armchair directions shows long 
range Neel correlations, and the lowest excitation spec- 
trum becomes gaplcss at the F point (for a finite size 
system there is a small gap which disappears after the 
extrapolation). For J2/J1 = 0.36, the SSCF in both 
of the zigzag and armchair directions shows short range 
Neel correlations, and the minimum of the lowest excita- 
tion spectrum remains at the F point, however, there is a 
large gap which does not disappear after the extrapola- 
tion. For J-2.1 J\ = 0.38, the SSCF does not show any long 
range correlation, and the short range correlations are 
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FIG. 9. (Color online) Momentum dependence of the lowest excitation spectrum for a system of size A'' = 2 x 50 x 50 
corresponding to the four different phases: (a)J2/Ji = 0.18 (Neel), (b)J2/Ji = 0.36 (GSL), (c)J2/Ji — 0.38 (staggered dimer 
VBC), and (d)J2/Ji = 0.48 (spiral). The dashed hexagon denotes the first Brillouin zone of the lattice. 



different along the zigzag or armchair directions, which 
is a indication that the lattice rotational symmetry is 
broken. Simultaneously, the minimum of the lowest exci- 
tation spectrum is away from the F point and the lattice 
rotational symmetry is clearly broken. There is also a 
gap in this region which remains finite in the thermody- 
namic limit. For J2/ Ji = 0.48, the SSCF shows magnetic 
long range correlations in both of the zigzag and armchair 
directions. Since one component of the ordering vector 
Qx = (corresponding to fc* = in the lowest excita- 
tion spectrum), the SSCF is Neel- like along the zigzag 
directions. This result agrees well with the spin wave 
calculations by Mulder et al.i^ 

Finally, we would like to talk about the next step of 
our work. We have used a mean field approach based 
in the Schwinger boson representation of the spin oper- 
ators. This mean field approach has the drawback of 
being defined in a constrained bosonic space, with un- 
physical configurations being allowed if this constraint 
is treated as an average restriction. This drawback can 
be in principle corrected by including local fluctuations 
of the bosonic chemical potential)^ This correction was 
calculated by Trumper et al^ for the Ji — J2 square lat- 
tice using collective coordinate methods, where a com- 



parison between the mean field results and the corrected 
results was made. However, this hard calculation allows 
only to calculate some special quantities like the ground 
state energy or spin stiffness. The corrections developed 
by Trumper et al. could be extended to spiral phases2£, 
which would allow to investigate, for instance, the present 
model. 



V. SUMMARY AND DISCUSSION 

In the present paper, we have investigated the quan- 
tum Ji — J2 Heisenberg model on the honeycomb lat- 
tice within a rotationally invariant version of SBMFT. In 
the region J2/J1 € [0,0.5], the quantum phase diagram 
of the model displays four different regionsj^ The mag- 
netic long range order of Neel and spiral types is found 
for J2/J1 < 0.2075 and J2/J1 > 0.398, respectively. For 
the spiral region, we get the spiral order from quantum 
disorder selection which agrees with Mulder et ali^ us- 
ing spin wave theory. In the intermediate region, the 
energy gap is finite while the local magnetization is zero, 
which indicates the presence of a magnetically disordered 
ground state. We have used the Z3 directional symmetry 
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breaking order parameter defined in Eq. (j32p to elas- 
sify this part into two different magnetically disordered 
phases: one is a GSL which shows short-range Neel cor- 
relations (J2/J1 S 0.3732), the other is staggered dimer 
VBC (lattice nematic), which breaks the Z3 directional 
symmetry (J2/J1 > 0.3732). Considering the properties 
of order parameters and the ground state energy, these 
three quantum phase transitions seem to be continuous. 

As we have mentioned above, recent theoretical studies 
of the phase diagram of the spin- 1/2 Ji — J2 Heiscnberg 
model on the honeycomb lattice have obtained conflicting 
results. The central controversial point is the existence 
and nature of magnetically disordered phases when the 
Neel order becomes unstable as increasing the frustration 
J2/J1. There is a growing consensus^2i33^i36^^^^^ 
that a magnetically disordered region should appear. 
However, the nature of this region is still not clear with 
different approaches giving different results. An early ED 
work by Fouet et alM^ flrst claimed that a GSL might ap- 
pear in the region J2/J1 ~ 0.3 — 0.35, and for J2/J1 ~ 0.4 
the system might be in favor of the staggered dimer VBC. 
A recent ED study by Mosadeq et alM- has claimed that 
a plaquette valence bond crystal (PVBC) might exist 
in the region 0.2 < J2/J1 < 0.3, and a phase transi- 
tion from PVBC to the staggered dimer VBC exists at 
a point of the region 0.35 < J2/J1 < 0.4. However, a 
more recent ED work by Albuquerque et al^, which has 
treated larger system sizes, has been unable to discrim- 
inate whether this magnetically disordered region cor- 
responds to PVBC with a small order parameter or a 
GSL. It is possible that the PVBC may just come from 
the finite size effectsi^ For larger J2/J1, it has been also 
hard to discriminate the staggered dimer VBC with spi- 
ral phases, since ED is especially difficult to treat the 
incommensurate behavior of spin correlations due to the 
small lattice sizes. 

There are two recent studies of this model using VMC 
with different variational wave functions. Clark et a/.— 
have used Huse-Elser states and resonating valence bond 
(RVB) states, and claimed that a GSL appears in the 
region 0.08 < J2/J1 < 0.3; a dimerized state which 
breaks lattice rotational symmetry for J2/ Ji si 0.3. How- 
ever, a more recent work by Mezzacapo et ali^ using an 
entangled-plaquette variational (EPV) ansatz have ob- 
tained lower energy estimates, and claimed that in the 
magnetically disordered region 0.2 < J2/J1 < 0.4, the 
PVBC order parameter vanishes in the thermodynamic 
limit. Therefore, the PVBC may just come from the 
finite size effects. Since the directional symmetry 
breaking order parameter has not been considered in this 
paper, it is still not clear that the lattice rotational sym- 



metry is broken or not in the region 0.2 < J2/J1 ^ 0.4. 

In a recent study using PFFRGi^ the authors have ob- 
tained that within the magnetically disordered region, for 
larger J2 / Ji , there is a strong tendency for the staggered 
dimer ordering; for low J2/ Ji , both of plaquette and stag- 
gered dimer responses are very weak. A further recent 
study using CCM^ has got a more quantitative magnet- 
ically disordered region: 0.207±0.003 < J2/J1 < 0.385± 
0.010, in which the PVBC phase has been reported. How- 
ever, the ground state within 0.21 < J2/ Ji ^ 0.24 is hard 
to be determined using this approach. 

The other controversial point is the form of the mag- 
netic long range order when J2/J1 exceeds the magnet- 
ically disordered region. There are two proposals: the 
anti-Ncel order— or the spiral order. It is difficult to get 
a conclusion by ED since it is hard to treat the incom- 
mensurate spin correlations due to small lattice sizes.— 
Both of the recent SE^ii and PFFRGi^ studies have not 
found any evidence for the existence of the anti-Neel or- 
der and concluded that the spiral state should be the 
stable ground state. However, both of the VMC with 
EPV ansata^^ and the CCM22 studies support the oppo- 
site proposal. Since we are interested in the exotic disor- 
dered phases in the magnetically disordered region and 
focus on J2/ Ji <= [0, 0.5], we can not exclude the possibil- 
ity that the anti-Neel order state exists for J2/J1 > 0.5. 

Due to the existence of strong quantum fluctuations 
and frustration, the spin-1/2 Ji — J2 Heisenberg model 
on the honeycomb lattice is a challenging model which 
needs further investigation especially for the nature of 
the intermediate phase. Unbiased numerical simulations 
are still needed, such as the density matrix renormaliza- 
tion group (DMRG) method,2i"I^ Recently, DMRG has 
been applied to spin-1/2 Kagome Heisenberg model^iiS 
and square Ji — J2 Heisenberg model2^, and obtained 
GSLs as the ground state. Since quantum fluctuations 
are expected to be stronger on the honeycomb lattice 
than those on the square lattice, it would be very inter- 
esting to apply DMRG to the spin-1/2 Ji — J2 Heisenberg 
model on the honeycomb lattice. 
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